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Abstract 


This thesis deals with the topic of generalized finite element method (GFEM) with particular 
reference to linear elastic fracture mechanics (LEFM). The presence of a crack in the domain 
leads to singular stresses at the crack tip, for linear elastic analysis. In this case, the conven- 
tional finite element analysis is known to lead to inferior quality crack parameters. The quality 
of extracted crack parameters can be improved only by significant mesh refinements at the crack 
tip. Use of refined meshes increases the computational load. In the present study, an effort 
has been made to characterise the crack parameters more accurately and (computationally) effi- 
ciently. The basic idea is to use specialized basis functions at the crack tip, using the asymptotic 
crack solutions. The method uses the partition of unity approach to make the approximation 
conforming. The code has been benchmarked with a large set of sample problems. Further, a 
crack propagation problem in planar elasticity is considered. Lastly, the problem of propagation 
of multiple, arbitrarily placed, cracks is considered. 
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Chapter 1 
Introduction 


1.1 Introduction 

The generalized finite element method (GFEM) is a combination of classical finite element 
method (FEM) and partition of unity method (PUM). In generalized finite element (GFEM), 
the standard finite element spaces are augmented by adding special functions which reflect the 
known information about the boundary value problem and input data; e.g. the singular solution 
obtained from the local asymptotic expansion of the exact solution in the neighbourhood of a 
centre point, etc. The essential idea in this method is to add enrichment functions to the finite 
element approximation which contains a discontinuous displacement field. The method exploits 
the partition of unity property of finite elements which was noted by Melenk and Babuska [1], 
and Duarte and Oden [2], namely that the sum of the shape functions must be unity. The special 
functions are multiplied with the partition of unity to construct an augmented conforming finite 
element space. In this way the local approximability afforded by the special functions is included 
in the standard finite element approximation, while maintaining the existing infrastructure of 
finite element codes. 

1.2 Literature survey 

In the field of computational mechanics, extensive research has been done in the effective use of 
the finite element method for problems in fracture mechanics. Effort has been primarily focussed 
on the modification of the approximation, incorporating the singularity present in the domain. 
Of the various elements developed, collapsed quadrilateral elements and quarter point element 
requires special mention [16]. While direct application of standard elements is the most straight- 
forward, a high degree of mesh refinement is required near the crack tip to capture the singular 
stress fields. 

Recently, various methods that attempt to do away with the mesh have become popular for 
solving boundary value problems, for example, the element-free galerkin methods of Belytschko 
et al. [11]; the hp Clouds method of Duarte and Oden [2]. These methods in their purest form 
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get rid of the finite element method (FEM) completely, and intend to replace traditional finite 
element codes with a new infrastructure. The methods have had limited success because they 
have not successfully addressed the problem of the numerical integration of the entries of the 
discrete equations. For example, all meshless methods must still employ an integration mesh for 
the numerical integration. However, because the regions of integration which are often intersect- 
ing disks in two-dimensions (2D) or spheres in three-dimensions (3D) do not coincide with the 
integration mesh, it is very difficult to control the error of the numerical integration. Further, the 
various meshless methods also require special techniques for applying essential boundary condi- 
tions, e.g., the method of Lagrange multipliers which raise the issues of stability. An enriched 
element free galerkin method (EFG) was proposed by Fleming et al [4]. However, it has yielded 
good results only when adequate refinement (in the background mesh) is used near the crack tip. 
Moreover, the stresses computed also tend to be oscillatory near the crack tip unless substantial 
refinement is used. All these modern meshless methods are based implicitly or explicitly on using 
a partition of unity over the the domain to ensure continuity of the approximation. The partition 
of unity method (PUM) was developed by Babuska et al [1]. 

Recently, the generalized finite element method (GFEM) was proposed by Strouboulis et al [6] 
to do away with the problems encountered in common meshless methods. The Computational 
Mechanics Company (COMCO) is currently developing a new fast and robust solver for struc- 
tural mechanics problem. It is based on a new solver technology, which promises to combine the 
modeling advantages and flexibility of finite element methods with the ease of use of meshless 
methods and with computational performance compatible or faster than modern finite element 
codes. The solver is based on proprietary extensions of meshless cloud methods and generalized 
finite element method (GFEM). It is being implemented in an hp-adaptive context within the 
computational environment PHLEX. Formally, the method is named generalized finite element 
partition of unity method or, for short, an element partition method (EPM) [7]. 

The present work, is a first step towards developing a generalized finite element based analysis 
platform for the study of quasi-static and dynamic crack propagation. The local asymptotic 
solution, at the crack tip, is chosen for enrichment. These special functions are made conforming 
by employing partition of unity based cut-off functions. 

1.3 Thesis organization 

The thesis is organised as follows: 

1. In chapter 2, we briefly review the concepts of the re-entrant corner and its generalisation 
to linear elastic fracture mechanics. The procedure employed to calculate the different 
crack parameters are outlined in this chapter. 

2. In chapter 3, a general procedure is implemented, based on energy principle, for deriving 
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the weak formulation studied in this work. 

3. In chapter 4, we discuss the main features of the generalized finite element method employed 
in this work. An adaptive integration scheme tailored to the present work has been detailed 
elaborately. 

4. In chapter 5, we present the numerical results obtained and a qualitative analysis has been 
done to compare the classical finite element solution and generalized finite element solution. 
Lastly an effort is made to predict the crack path using the above principle. 

5. Finally in chapter 6, we summarize our work and draw the relevant conclusions and the 
scope for future work. 
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Chapter 2 


Linear elastic fracture mechanics 

2.1 Two dimensional planar elasticity 

Equilibrium Equation 

In planar elasticity problems, the variation of stress components in third direction is negligible. 
Also the body force part in fracture mechanics is absent, hence the equilibrium equations can be 
simplified as follows 


d<Jn dai2 
dxi dx2 

dcr 12 dcf22 

dxi dx2 


( 2 . 1 ) 

( 2 . 2 ) 


Strain displacement relationship 


If Ui and U 2 are the two planar displacement components, then the strain-displacement relation 

is given by 


dui 

£^11 — a 

UXl 

(2.3) 

dU2 
- dx. 

(2.4) 

1 / dui du2 N 

(2.5) 

2\dx2'^ dxi) 

The compatibility relation is given by 


d^su d‘^S22 ^ d‘^ei2 _ 

dxi dx\ ^dxidx 2 

( 2 . 6 ) 


Stress strain relationship 

For the linear elasticity problem with isotropic materials, the stress strain relationship can be 
given as 

£n = - ^{<^22 + ^ 33 )] (2.7) 
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Plane deformation 


£22 = 

^[0-22 — + eras)] 

(2.8) 

£^33 = 

^[<^33 ~ i^(c’'ii + £^ 22 )] 

(2.9) 

£■12 = 

(l + v) 1 

E " 2G 

(2.10) 


For the case of plane stress i.e. for a very thin plate, the out of plane stresses are negligible in 

the interior points of the plate. Thus we have 


(7i3 = (723 = <733 = 0 

( 2 . 11 ) 

Thus the stress strain relation reduces to 


£11 = - 2 ^ 0 - 22 ] 

( 2 . 12 ) 

£22 = '^[‘^22 — iso’ll] 

(2.13) 

1 

ei 2 - ^^12 

(2.14) 

Similarly for the plane strain case, i.e. for sufficiently thick plates 
strain leads to 

, the through the thickness 

£13 = £23 = £33 = 0 

(2.16) 

c’33 = + C 22 ) 

(2.16) 

Thus the stress strain relation reduces to 


(1 - 1^2) r u ^ 

£11 - g ku 

(2.17) 

(1 - 1/2) r u 

£22 = g k22 

(2.18) 

(1 — 1/2) I- V , 

(2.19) 

or, a general representation for the planar constitutive relationship 

can be given by 

£11 = ^[<^11 — ^^*<^ 22 ] 

(2.20) 

£22 = ^[<722 - 

(2.21) 

1 + U* 

£12 - ^^12 

(2.22) 

where 


( E for plane stress 


E* = < 

(2.23) 

1 for plane strain 
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V 


• V for planestress 


(2.24) 


( pip for planestrain 

2.2 Linear elastic fracture mechanics 

Below, a brief outline of the theory of fracture mechanics for linear, plane elasticity, is given 
Consider a vertex A, as shown in figure ( 2 . 1 ). In the case of planar elasticity, the Airy stress 



Figure 2 . 1 : An arbitrary domain with a re entrant corner. The local cartesian and polar coordi- 
nate system with respect to the vertex A 


function U satisfies the biharmonic equation : 

/ 5^ 15 i 5^ \ /_5^ 1 1 r r _ Q 

v5r^ r 5r 59 ^) v5r^ r 5r 89 '^) 

Let us consider solutions of the form (see [ 12 ]) : 

U = U{r,9) = r^+^F{9) 


(2.25) 


(2.26) 


Therefore , F{9) must satisfy : 


F"" 2(A2 -1- 1)F" + (A^ - 1)F = 0 


(2.27) 


where the prime represent differentiation with respect to 9. The general solution of equa- 
tion (2,27) for A ^ 0 and A ^ ±1 is : 

F(0) = Oi cos(A — 1)^ -h 02 cos(A -I- 1)^ + 03 sin(A — 1)0 H- 04 sin(A + 1)0 (2.28) 
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The modifications necessary for the cases A = 0 and A = ±1 are obvious. The stress components 
are related to the stress function U as follows : 


0~ ^ —— 
aQ = 

'^rO = 


1 dU 1 d'^U \ 
+ — — = r^' 


r dr 
d^U 

d^i 

1 d'^U 


r2 


'■[{X + l)F{6) + F"m 


= A(A + l)r 

2 de 


A- 


^F{e) 


r drdd r 


= -\r^-^F'{e) 


(2.29) 

(2.30) 

(2.31) 


Let us assume that the edges at ^ = ±Y are stress free which implies that : oq = = 0 

on these edges. Consequently, from equation (2.28), (2.31) and (2.31), after straightforward 
algebraic manipulation, we have 


cos(A — 

1)y cos(A-l-l)^ 

0 

0 


\ 

ai 

Asin(A - 

-l)f sin(A + l)f 

0 

0 

< 

^2 

> 

0 

0 

sin(A-l)f 

sin(A + l)y 


03 

0 

0 

AcOs(A — l)y COS(A + l)y_ 


. ^4 , 


(2.32) 


where A is defined as 


A - 1 
A + 


(2.33) 


This gives rise to two sets of decoupled simultaneous equations in terms of Oi, aa and as, a^. 
A non trivial solution exists only if the corresponding determinant vanishes. This is possible if 
either 


cos(A - 1)^ sin(A + 1)^ - Asin(A - 1)^ cos(A + 1)^ = 0 (2.34) 


or 

sin(A - l)^cos(A + 1)^ - Acos(A - l)^sin(A + 1)^ = 0 (2.35) 

To find non trivial solution for equation (2.32) we set as = 04 = 0 and find A such that equa- 
tion (2.34) satisfied, or we set ai = as = 0 and find A such that equation (2.35) is satisfied. 
Equations (2.34) and (2.35) can be simplified to 


sin Aa + A sin a = 0 (2.36) 

sin Aa — A sin a = 0 (2.37) 


Now A can be complex and either a simple or multiple root of equation (2.36) or equation (2.37) 
If A is complex, then its conjugate is also a root of equation (2.36) or equation (2.37) , and both 
the real and the imaginary parts of r^+^T’(0) are solutions of equation (2.25). If A is a solution, 
then -A is also a solution; however the corresponding stress field has finite strain energy only if 
the real part of A is greater than zero. 
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Assume that ^ (i— 1,2....) is a solution of equation (2.36) and is real and simple. Then 
from equation (2.32) : 


ai cos(Ap^ _ 1)^ + a 2 cos(Af ^ + 1)^ = 0 

sin(Ap^ ~ 1)^ + a2sin(Af^ + 1)^ = 0 
z z 

Let us define : 

gW cosjXP - l)f ^ Af^ sin(Af^ - l)f 
cos(Af^ + l)f sin(Ap^ + l)f 

With this notation, equation (2.26) can be written as : 

u = r^‘'^+^(cos(Af ^ - 1)0 + gf ^ cos(Af ^ + 1)0) 


(2.38) 

(2.39) 

(2.40) 


(2.41) 


Similarly, if A^^^ is a real, simple root of equation (2.37), then the stress function is : 

U = r^^''^+^(sin(Af ^ - 1)0 + gf ^ sin(Ap^ + 1)0) 

where Qf^ is defined as : 

( 2 ) ^ sin(Af^ - l)f ^ Ap^ - 1 cos(aP - l)f 
sin(Af^ + l)^ Ap^ + 1 cos(Af^ + 1 )y 

Equation (2.41) is symmetric with respect to 0 and equation (2.42) is anti symmetric with 
respect to 0. From equation (2.41) and equation (2.42), expressions for stress and strain can be 
derived. The stress and displacement field corresponding to equation (2.41) are called Mode 1 
fields which are given by : 


(2.42) 

(2.43) 


Uii = '[(« - (5P^(At(l) + 1)) cos Ap^0 - Ap^ cos(Ap^ - 2)0] 

« 2 p = '[(« + <5P^(Ai(l) + 1)) sin Ap^0 + Ap^ sin(Ap^ - 2)0] 

which can be written as 



where k is Kolosov’s constant given by ; 


(2.44) 

(2.45) 


(2.46) 


K = 3 — 4i/ for plane strain 

K = - — - for plane stress 
l + iy 

The Mode I stress tensor components are ; 


ap) = ApV^>-'’-p2 - gp) (aP^ + 1) cos(Ap^ - 1)0 - (Ap^ - 1) cos(Ap) - 3)0] 
= ApV^‘-'’"^[2 + QPpAp^ + 1) cos(Ap^ - 1)0 + (Ap^ - 1) cos(Ap^ - 3)0] 
= ApV^*' ^■p(Ap^ - 1) sin(Ap^ - 3)0 + Qp^(Ap^ + 1) sin(Ap^ - 1)0] 


(2.47) 

(2.48) 


(2.49) 

(2.50) 

(2.51) 
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The stress and displacement fields corresponding to equation (2.42) are called Mode II fields 
which are given by : 

4? = - QfHAi(2) + 1)) sin - Af’ sin(Af^ - 2)6] (2.52) 

4? = + Qf^W) + 1)) cos \f^e + Af’ cos(Af^ - 2)9] (2.53) 

which can be written as 

uf’ s I JJ I = (2.54) 

The Mode II stress tensor components are : 

= AfV^®-i[2-QP(Af) + l)sin(Af -l)0-(Af)-l)sin(Af)-3)0] (2.55) 

42 = [2 + Qf^ (Af ) + 1) sin(Af ) - 1)0 + (Af ^ - 1) sin(AS'’ - 3)0] (2.56) 

<7i 2 = -Af V^'''-'[(Af - 1) cos(Af - 3)0 + Qf (Af ^ + 1) cos(Af ^ - 1)0] (2.57) 


Thus in the general case in the vicinity of the re-entrant corner, the solution can be written using 
equation (2.46) and (2.54) as : 


{“> = E + E 


2G ^ * V-" ' ^ 2G 

i=l i=l 


(2.58) 


The coeflficients ^4^ ^ and Af^ are called generalized stress intensity factors. In case of a crack we 
have a = 27r. Equations (2.36) and (2.37) are identical (sin27rA = 0), and all roots are real 
and simple : 


A 


( 1 ) 


Ap^ = ±7r, ±n, ±2, ±- 


(2.59) 


The coefficients ^ and Af'^ are related to Mode I and Mode II stress intensity factors of 
linear elastic fracture mechanics. In linear elastic fracture mechanics three modes of crack opening 
are present in three dimension. In plane elasticity, only first two modes are present. These are : 


1. Mode I - Tensile Mode or Opening Mode 


2. Mode II - Shear Mode or Sliding Mode 

Below, we give the expression for the displacements and stresses for the mentioned two modes. 


Mode I 

For Mode I type of loading 


the solution in the vicinity of the crack tip is given by 
Ki 

- r.os . 

2 ‘ 


0-11 


cos -[1 


a'22 — 


a'12 — 


. 0 . 30 , 

, ^ _ sin - sm 

•s/2'Kr 2 2 2 

Ki 0 ., . 0 . 30 , 

^ccSjU + smjsm-l 

Ki 9 . 6 39 

2 2 2 


(2.60) 

(2.61) 

(2.62) 
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where Kj is the stress intensity factor in Mode I. The displacement components Ui and U 2 are 
given as follows 


Ui 

U2 


Kj f r \ k Or . , 

■gIsJ 2 [1 - 2" + sm- - 

1 


f r \ h 9 \ 

= 'gIsJ sm-[2-2i. + cos 


e_- 

2 . 


(2.63) 

(2.64) 


(733 = 0 for plane stress 

(^33 = ^'(<711+ < 722 ) for plane strain 


Mode II 


For mode II type of loading the solution in the vicinity of crack tip is given by 


Kij .9 9 30 

022 = — F= sm - cos - cos — 
^/^ 2 2 2 

Kii 9 . 0 . 30, 


(2.65) 

( 2 . 66 ) 
(2.67) 


where Ku is the stress intensity factor in Mode II. The displacement components Ui and U 2 
are given as follows 


Ui 


9 


_ Ku ( r \ \ ^ 

G V27r/ ^'''21 


2 — 2z/ + cos^ - 


Kji ( r \ 2 9 

U 2 = I — 1 cos X 


-1 + 2n + sin^ x 


( 2 . 68 ) 

As it is evident from the functional form of (7^ that singularities in the form of are present. 
This means that when r — )■ 0, an and 0-22 becomes infinite. 


Mixed mode analysis 

The near tip displacement fields for combined Mode I and Mode II loading are obtained by 
superimposing the two fields and is given by : 


Ki ( r \ \ 0, ^• 2^1 { T \k . 9 . , „ 2^1 

— ) cos-[«;-l + 2sm -] + -^ ( — ) sin -[« + 1 + 2cos2 -] 


U 2 


2G\2ttJ ’"2 
Ki fr\\ . 9 
2G \2 tt) ®^^2 

Similarly, the stress field is given by ; 

^ cosf( 

0 


0 , 


2G \ 27 rJ 
Kn / r \ h 


+ 2GV2^ 


2%) 


4 0-2^3 

COS -\k — 1 — 2 sin - 
2^ 2^ 



Kr 


1 — sin I sin ^ 


\/27rr 


cos f + sin X sin 


0.;. 30' 


Kn 


cos f sin ^ cos ^ 




7rr 


— sin ^(2 + cos I cos 
sin I cos I cos ~ 


f) 

cos |- ( 1 — sin sin 


0 


inf) 


(2.70) 

(2.71) 


(2.72) 
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The above general definition of the asymptotic solution near the crack tip is fundamental to 
the process of quasi static crack propagation. Below, we define the J integral that is used to 
characterise the criticallity of an existing crack. 


2.2.1 J integral 


Taken along a material curve F which is traversed in the counter clockwise direction between the 
two crack sides , the J integral is defined as 


J = 



(2.73) 


where W = aijdsij is the deformation energy per unit volume, s is the arc length and Ti is 
the traction exerted on the body bounded by the curve F and crack surface. 


T) — O'jjTZj 


where rij is the unit outward normal to the curve F. Two important characteristics of J integral 
^2 



Figure 2.2: Computation of J integral 


are 

• It is path independent for linearly elastic bodies. 

• For lin early elastic bodies, J integral represents the energy release rate and is same as the 
energy release rate G, proposed by Griffith [15]. 


12 



It can be shown for Mode I loading condition 


and for Mode II loading condition 
Thus for mixed mode analysis, 


J 


J = 


E* 

E* 


Kr^ , Kf 


J = ^ ^ (2.74) 

E* E* 

Thus we note that the computation of the stress intensity factors is essential for obtaining the J 
integral 

Computation of stress intensity factors 

The stress intensity factor can be computed by several methods. Below we consider one of 
the method outlined in [3], which employ contour form of the interaction integral. The local 
coordinates are the crack tip coordinates with the Xi axis parallel to the crack faces (see figure 
2.2). Two states of cracked body are considered following [3]. State 1 corresponds 

to the actual state and state 2 {alpslpuf) is an auxiliary state which will be chosen as tlie 
asymptotic field for Mode I and Mode II. The J integral for the sum of the two states is 

d{u}+ujy 


j(i+2) ^ ^ - (al + afj) 

Expanding and rearranging the terms give 

j(l+2) ^ j(l) ^ j(2) ^ j{l,2) 


rijdT 


where is the interaction integral for states 1 and 2, 


/(b2) = - al 

and W^h2) jg interaction energy, given by 


(1) duf (2) du] 1 


— a} 


dx, dx.y^^ 


riidV 




Substituting equation (2.74) for the combined states and rearranging the term gives 


j(n-2) = J(1) + J(2) + + K\yK\f) 

Equating equation (2.76) and (2.79) leads to the following eqaution : 


(2.75) 


(2.76) 


(2.77) 


(2.78) 


(2.79) 


(2.80) 
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Now choosing state 2 as Mode I asymptotic field with Kf'^ = 1 and = 0 gives Mode I 
stress intensity factor for state 1 or actual state in terms of interaction integral as 

I) ^ 2 . 81 ) 

2 

Then choosing the auxiliary state as that for Mode II asymptotic field with Kf'"^ = 0 and 
~ ^ gives Mode II stress intensity factor for state 1 in terms of interaction integral as 

Kf} = EljiiMode II) ^ 2 . 82 ) 

2.2.2 Crack growth criteria 

We briefly review the crack growth law used to specify the direction of crack growth. The criteria 
for determining the crack growth direction can be summed up in the following three categories : 

1. The maximum energy release rate criterion 

2. The maximum tangential stress criterion 

3. The minimum strain energy criterion 


In the present work, we used the maximum tangential stress criterion or the maximum circum- 
ferential (hoop) stress criterion, which is identical to maximum energy release rate criterion for 
linear elastic fracture mechanics. The maximum tangential stress criterion states that the crack 
will propagate from its tip in a direction 6c so that the circumferential stress agg is maximum. 
Therefore, the critical angle 6c deflning the radial direction of propagation can be determined by 
setting the shear stresses associated to zero because this corresponds to the principal direction. 
Thus 




^ '^^Ki sin 6 + \:Kii (3 cos 6 — 1) 


(2.83) 


This leads to the equation defining the angle of crack propagation (in the tip coordinate system) 


6c 


Ki sin 6 c + i^//(3 cos 0c — 1) = 0 


(2.84) 


Solving the above equation gives 

e, = 2arctani(^±y(^f + 8) 


(2.85) 


For Kii = 0, 6c = 0, i.e. the crack propagates straight, without turning. 
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Chapter 3 

Finite element formulation 


3.1 Finite element formulation using energy principle 

The total potential for the plate is given by 

N 

n (u) = (u), (3.1) 


e=l 


where tt® is the total potential of the non-intersecting (but adjacent) sub-domains e which are 
part of the domain (N sub-domains are considered here). The total potential can be expressed 
in terms of internal strain energy and external work done as follows: 

TT® (u) = [/(^) - TF^®^ (3.2) 

n(u) = ^ / Wxx{u)SxxivL)+axy{n)'rxyi^)+(^yy{-o.)€yy{v.)}dA- Jj^f^u+fyv)dA- J^{TxU+Tyv)ds 

(3.3) 

where {f} is the body force vector, {T} is the traction vector and {u} is the displacement vector 
given by : 

« = { t } 

= { ?: } 

Minimizing the total potential energy gives 

n( u-|- gw) - n(u) _ Q 


d^^)(n(u)) = lima-^Q- 


a 


where {w} is the perturbation vector 


r 'i / '^X 

{-} = I 
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5(^^(n(u)) - / (cr2:i(u)£:2;a;(w) + crj.y(u)7(w) + £r 2 /y(u)£j/ 3 /(w)) - f (f^Wx + fyWy)dA 


{TxWx + Ty'Wy)dS 


Hence the bilinear form can be written as : 

/ 3 (u,w) = icrxx{u.)exx{'v^) + crxy{u)jxyi'^) + ayy{u)eyy{-w))dA 

J A 

■^(w) = J {fxVJx + fyWy)dA + j^{TxWx+TyWy)dS 

3.2 Finite element approximation 


(3.4) 

(3.5) 






ELi Ui^i{x,y) 

YA=iVi^iix,y) 


0 $2 3 ^3 0 

0 $1 0 $2 0 $3 


Ui 

Vi 

U2 

V2 

U3 

V3 




Similarly, 


{w} 

Ei=l 

{w} = mw} 

Triangular elements are used in the finite element approximation employed in this study, along 
with hierarchic shape functions of order p (p < 3). rz is the number of degrees of freedom for 
each element and it is equal to for triangular elements. are shape functions in global 

coordinate system x and y. 

Let the stress and strain vectors corresponding to {u} be {cr} and {e} which can be defined as: 


M = <j a, 

'^xy 


{4 = 


(3.6) 


Jxy 


Prom the generalized Hooke’s law which relates stress components to the respective strain com- 
ponents in global coordinate system, we have 

W = [Q1{4 (3.7) 


1R 



where the material stiffness matrix [Q] is given by : 

^ I'S n 

1—1/2 '•J 

1QI= rSs 0 

.0 0 G 

Ejl-v) vE 


for plane stress 


1QJ = 


Ed— I 


uE 


for plane strain 


0 0 c? _ 

The strain can be represented as: 

{e} = [D]{u} 

where [D] is a differential operator in terms of global coordinates such that 


[D] = 


^ -k 0 


0 # 

dy 


JL A 

dy dx J 


Hence, the strain field using finite element approximation is given by ; 

= [DJIujTiJAr} 

Substituting the expression for mfem, we have 

{e(Uir£;Af)} = [B]{U} 

where [B] is given as : 


b; = 


Similarly, for 

£(w) = [B]{W} 

Thus, in finite element approximation the bilinear form can be written as : 

= j {e’’ [v)){cr(UFEM)]iA 
/5(ueem,w) = |B]{W})^[D]([B]{U} 


■ d^ 

0 


0 

d^ 

0 .. 

- 

dx 

dx 

dx 


0 


0 

d^9. 

0 

3^ 


dy 

dy 

dy 




a$2 



d^9 


- dy 

dx 

dy 

dx 

dy 

dx 

*' » 


(3.8) 


(3.9) 


(3.10) 


(3.11) 


(3.12) 


(3.13) 


(3.14) 
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(3.15) 


w) = {W}^[B]^[D][B]{U} 
The linear functional can be written as : 


F{w) = j fwdA + J^^Twds 


F{w) = {W}^ [ J + J_^[^f{T}dsj (3.16) 

Now, substituting the above expressions for bilinear and linear functional , w^e get 

= F{w) 

{W}^[^[B]^[D][B]d^]{U} = {Wf[J[^f{f}dA + J^[^f{T}ds 

Since, {W} is arbitrary, we have 

[K]{U} = {F} (3.17) 

where [K] is the element stiffness matrix given by ; 

[K] = f [Bl^|D)lB]dyl (3.18) 

Ja 

and {F} is the load vector given by : 

{F} = [ [#]^{f}ds4 + [ [<&]^{T}ds (3.19) 

J A JV 

3.2.1 Computation of stiffness matrix and load vector 

After expanding the terms of equation (3.18), the entries for element stiffness matrix can be 
written in the following form for a given i and j : 


K(2i-l){2j-l) 

K(^2i-l)(2j) 

A'(2i)(2i-1) 

K(2i){2j) 


/ (Oil + Q33^j,y^i,y)dA 

JA 

j (0l2^j,j/’^i,a; T Q33^j,x^i,y)dA 

JA 

J" {,Q\2^ j,x^i,y T Q33^ j.,y^i,x)dA 
[ iQ22^j,y^i,y + Q33^j,x^i,x)dA 


Here z and j varies from 1 to iVT)OF (where NDOF is the number of degrees of freedom). 
Similarly, from equation (3.19) we get 

F(2i—1) ~ J fx^idA + Tx^ids 


F< 


(2i) 


fy^idA + 


Ty^ids 
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3.3 Geometry mapping 

Numerical integration of the stiffness matrix entries and the load vector entries necessitates the 
transformation of the physical coordinates (x, y) to the master element coordinates (^, ??). Hence 
instead of defining the shape functions in the physical coordinate system, it is defined in master 
element coordinates denoted by Nj. The geometry is expressed in terms of shape functions : 

3 

(3.20) 

2=1 
3 

i=l 

where xf and yf are the physical coordinates of the i-th node of the element e and Ni is the shape 
functions defined in master element coordinates (^,7?). 

Linear mapping 


X = 

y = 




Figure 3.1: Linear mapping 

Linear mapping is generally used when the sides of the physical elements are straight edges. 
Then the shape functions are given as : 

Ni = 1-i-V 

n 2 = e 

N3 = V 

Substituting the above relations, we get 

X = xl{l - ^ - rj) + xl^ + xIt} 

y = ytii - ^ - n) + y2( + yzv 

or 
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( 3 . 22 ) 


{ X — x\ 

I y-vl 


{xi - x\) {xi - x®) 1 r ^ 1 
(yl - y{) (yt -yt) \ \ y j 


/ ^ \ ^ f iyl-yf) -(4 - xf) 
Ivj |J|L-(y|-yO 

where [J] is the Jacobian matrix given by : 


(x-xf) 1 

(y-yf) J 


[j] = 


(xl-xf) (xl-xf) 

(yl - yt) (yl - y!) 


The derivatives of the shape functions are calculated using chain rule : 


d^i _ dNjdT] 

dx dx dr] dx 

d^i _ dNjdr] 

dy d^ dy dr] dy 


The derivatives for linear transformation are given by : 


dx 

dy 

dr] 

dx 

dr] 

dy 




( 3 . 23 ) 


( 3 . 24 ) 

( 3 . 25 ) 

( 3 . 26 ) 

( 3 . 27 ) 


Numerical integration 

The numerical integration over a general triangle is achieved by mapping the triangle to a defined 
master element and then performing the area integral. The mapping can be either linear or 
quadratic or any order which depends on the geometry of the general element is described in the 
above paragraph. After transformation, the integral value is given by : 

[ F{x,y)dxdy = [ F{i,r])\J\d^dr] 

Ja^ Ja 

. MINT 

- 2 E 

i=l 

where A® is the area of the physical element, A is the area of the master element and NINT is 
the number of integration points. 
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Chapter 4 

Generalized finite element formulation 


4.1 Introduction 

In the standard finite element method, the approximate solution is constructed by employing 
piecewise polynomial function of degree p. A major drawback of this method is that for the 
domain with re-entrant corner (or crack), severe mesh refinement, near the vertex of the corner 
is required. Thus the size of computational problem increases. This is because the solution is 
very unsmooth in the vicinity of the crack tip. However in case of linearised elasticity based 
fracture mechanics, the representation of the unsmooth part of the solution is known in the 
vicinity of the crack tip, as a series. Hence, incorporating this special function in the finite 
element representation will significantly improve the quality of the finite element solution in 
the vicinity of the crack tip. Thus the necessity of the mesh refinement would be drastically 
reduced for this type of finite element analysis. Further, the critical parameter (for example the 
J integral) can be obtained more accuartely as compared to conventional finite element analysis. 
Below, a detailed outline of the generalized finite element method (GFEM) is given. 

4.2 Construction of the conforming generalized finite el- 
ement approximation 

In this section, we describe how the finite element approximation can be generalized to include 
special functions. The standard shape functions used in finite elements result in mapped poly- 
nomials and are known to approximate smooth functions well. However, if the solution is known 
to be harmonic, it would be advantageous to use harmonic polynomials to form the basis. But 
the problem of continuity of approximation or conformity should be addressed properly to make 
sure that the approximation has finite energy. Conformity has been addressed in the standard 
finite element method (FEM) by constructing element shape functions in such a way that the 
proper number of degrees of freedom match at the boundaries of neighbouring elements and thus 
create a continuous approximation. Using any of the standard FEM elements and an associated 
family of shape functions of degree p, we can exactly represent all polynomials of degree p. These 
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shape functions are pasted together from neighbouring elements and form a basis by suppressing 
any linear dependencies, e.g., rigid body modes in elasticity. The resulting system of equations 
are solvable. The partition of unity method allows creation of conforming approximation using 



Figure 4.1; Linear support hat function 



Figure 4.2: Special function to be used for enrichment for e.g., 



Figure 4.3: Generalized basis function given by - 2x) over the domain (0, 0.5) 

special functions. Given a set of overlapping subdomains or patches {QJ (see figure 4.1) and a 
set of functions with desirable approximation properties associated with each patch = {'0}}) 
(see figure 4.2) we define the PUM approximate solution over the domain Q (see figure 4.3) as 

Upi/M = ^2 

i j 


22 



where ij)^j[x,y) is a special function from the set {■0*}. The space Tj is called the patch space. 
The a] are the constant coefficients. The sequence of functions (one for each patch of the 
set is a partition of unity on SI and serves to enforce continuity. 

Although the partition can be constructed in any manner as long as it satisfies the definition, 
in the present work we construct it using linear hat functions over coarse patches of triangular 
elements (see for one dimension example figure 4.1 - 4.3). The patches {Slj} and each associated 
member of the sequence {$j}must satisfy the conditions stated below : 


is an open cover of Cl, } 


e V i, 

(4.2) 

sup$j G closure{Cli) V i, 

(4.3) 

= 1 on Q, 

(4.4) 

m\T s c-oo, 

(4.5) 


The displacement field, at a crack tip, is contained in the span of the following four functions : 

( 9 6 6 6 ') 

{ 7 (r, 0)} = I \/r cos -, x/rsin -, -v/r sin - sin^, \/r cos - sin^ i (4.6) 


These functions will be used to enrich the standard finite element shape functions to give better 
approximation. Although more general form of enrichment function can be used, we used the 
above span to enrich the shape function used by Belytschko et al. [4]. Thus we can completely 
generalize the finite element approximation, allowing any family of special functions to be mixed 
with the mapped polynomial shape functions. The generalized finite element approximation is 
given by 

n-uert npEM 

'^GFEM = 

i=l j k 

where n„ert is the number of hat functions at element vertices, rii is the number of special functions 
in the patch space associated with the vertex and ufem is the number of finite element side 
modes and interior modes. 

Note that the traditional FEM basis functions are evaluated on the master element (^, t]) but the 
higher order special functions are expressed in terms of unmapped physical coordinate system of 
the problem {x,y). 


4.3 Linear dependency 

As mentioned earlier, the GFEM may lead to a singular stiffness matrix because of linear depen- 
dence of the local approximation functions. On the contrary, in traditional FEM the solution 
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corresponding to zero eigen value is removed by imposing essential boundary conditions. However 
this is not possible in GFEM because the eigen functions associated with the zero eigen value are 
not known generally. To solve such a linear system of equations one can use several methods like 
the perturbation approach [6]. In the present study, perturbation technique is incorporated in 
the solver code to solve such system of equations, though singularity does not arise in our case. 
When [A] is a singular matrix, then a particular solution can be obtained by : 

[^]a = [^] + a [I] 

where e is the perturbation, usually taken very small and [I] is the identity matrix. Thus [A]; 
is positive definite matrix and can be inverted by any standard method. 


4.4 Adaptive integration 


Since in the present work we are enriching the usual finite element space with local asymptotic 
solution in the vicinity of crack tip which have singular derivatives at r = 0, we must take care of 
numerically integrating their entries in element stiffness matrix or load vector. This is done using 
adaptive gauss quadrature integration, in the elements where the special functions are employed. 
For all other elements the standard gauss quadrature rule is applied. In the following paragraph, 
adaptive integration scheme will be discussed elaborately. 

The GFEM approximation will be defined as 

1 9 

= cos 

where Nmi^, v) is the traditional finite element shape function defined in the master element and 
Ni{^,ri) is the new auxiliary shape function. Now to compute the stiffness entries, we need to 
' find the derivatives of the auxiliary shape functions which we can find by applying chain rule. 


dNjj^, v) ^ ^ dNra{^,v) dC dNmi^, ri) dr] 
dx I dx dr] 9xJ 


r 2 cos 2 + 


. d I 9 


So, the stiffness entries corresponding to the auxiliary degree of freedom involves singular deriva- 
tives. To compute the stiffness entries for such row and column which involves the auxiliary 
degree of freedom, we divide the master element into four equal parts and compute the area 
integral on each of the sub master element and added to give the total value for that element 
(see figure 4.4). This value is then checked with the value calculated in the previous iteration, if 
the error is within the range of the tolerance specified in the code then no further subdivisions 
are carried out and calculation of next stiffness entry is carried out. However, if the error is more 
than the tolerance specified then out of the four sub master element, the sub master element 
asociated with the crack tip is redivided into four equal triangles and the integrand is evaluated 
on each of the element parts and summed to give the total value (see figure 4.4). Special care 
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should be taken to ensure that the weighting factor of the standard gauss quadrature rule gets 
modified after a triangle is divided. For example, in the first iteration after division of the master 
element into four equal triangles the weighting function for all the sub master elements should 
be divided by four. While the situation in any other iteration greater than one is more complex, 
since it involves both the sub master elements obtained in the previous iteration and the integral 
domain obtained after the division of the triangle associated with the crack tip into four more 
parts. Hence special care should be taken so that all the previously divided elements which have 
not been subdivided should retain their weighting function pertaining to the previous iteration 
and for the triangle which has been sub divided, its weighting function should be divided by four 
again to get the weighting function for each of the sub divided triangles. In the present work a 


3 Enriched Node 



Figure 4.4; Scheme for Adaptive Integration 

ninth order integration rule is applied which results in nineteen integration points. In each of the 
sub divided triangles, the same order integration rule is applied to maintain the simplicity of the 
code. An important aspect of the present adaptive integration procedure is that we apply the 
prior knowledge of the node corresponding to crack tip which has been enriched and subsequently 


25 



attack the node to achieve desired accuracy. This approach ensures that the stiffness entries are 
accurate. However one has to be cautious since the same integration rule is applied to each of the 
sub divided domain, proper care should be taken to transform the sub master element coordinate 
to the actual master element. This is due to the fact that the shape functions are defined in 
actual master element and not in the sub master element. This is taken care of by a linear map 
of the sub master element coordinates to the actual master element coordinates since the master 
element is an isoceles right angled traiangle having straight edges. 

An important observation is that the present approach takes longer time for diagonal elements. 
This is due to the fact that the off-diagonal elements exhibit weaker singularity than the diagonal 
elements. In this work, we imposed a tolerance of 10“°®. Note that, we could have carried out a 
single computation of unsmooth integrand, and obtain the desired level of subdivision required 
in the master element, for elements at the crack tip. This can then be used to do the integration 
in one shot, for all the elements at the crack tip. 

4.5 Solution technique 

The solution phase of the finite element computer programs consists of three parts: 

1. Assembly of the stiffness matrix and load vector from the element level stiffness matrices 
and load vectors; 

2. Enforcement of the principal boundary conditions, and 

. 3. Solution of the system of simultaneous equation. 

In engineering problems the size of the system of equations is generally quite large, and the entire 
global stiffness matrix cannot be stored in runtime memory, even if only the non-zero entries are 
stored. It is important to design the solution algorithm so that the arithmetic operations in the 
solution phase can be performed as efficiently as possible. 

In the present work, we have used the frontal solution method. The frontal solution method was 
proposed by Irons [12]. Its main feature is that the order in which the gaussian elemination is 
performed is element oriented rather than node oriented. It is not necessary to number nodes for 
minimum bandwidth. Another important feature of this method is that assembly of the stiffness 
matrix and the gaussian elimination is carried out concurrently. 
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Chapter 5 

Numerical results for stationary crack 


In the present chapter we will discuss the results obtained. Both straight and arbitrary shaped 
crack segments were considered. Initially a comparison is drawn between FEM solution and 
GFEM solution to prove the efl&ciency of the proposed methodology. The cases were solved using 
hierarchic shape functions of order p (p<3). This could be very easily extended to even higher 
order approximating polynomial. In all the cases, the enrichment functions used were conforming 
to the one proposed by Fleming et al [4]. However more general enrichment functions could be 
employed. 

Crack modelling 



Actual Crack 



Modelled Crack 


Figure 5.1: Crack modelling 

In the present work, any arbitrary shaped crack is modelled as a series of straight line segment 
connected at nodes as shown in figure 5.1. To ensure the discontinuity along the crack faces, two 
nodes were created at the same position one on the upper surface of the crack and the another 
one on the lower surface. This ensures that the element discretization above and below the crack 
faces are independent of each other. The local coordinates at the crack tip are taken parallel to 
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the end line segment (see figure 5.1). 

In all the cases the full plate analysis is done to elaborate the computational strategy in case of 
multiple singularities in the domain. For all the problems, plane strain condition is chosen and 
the following material properties, boundary conditions and loads are taken : 

Material properties : E = 100 kpsi and u = 0.3 
Boundary conditions : Four corner nodes fixed 
Traction ■. Ty = 1.0 psi on the opposite edges 



Figure 5.2: Boundary conditions and loading conditions 

The plate is taken to be of dimensions 5in x 3.5in (as shown in figure 5.2). 

Problem 1 : Centrally located straight edge small crack 
A center cracked plate with the following parameters is analysed : 

Domain = 5 x 3.5 ir? 

Crack Length = 0.5 in 

Two different mesh is employed, one is coarser(see figure 5.3) and the other one is finer(see figure 
5.4). The computed values of the various crack parameters using either the conventional FEM 
or GFEM are reported in table 5.1 and 5.2. Note that in the computation of stress intensity 
factors, the contour is taken to correspond to the outer boundaries of either one layer of elements 
adjacent to the crack tip or two layer of elements. 
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Figure 5.3: Straight edge crack of length 0.5in with coarse mesh 
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Table 5.1: Stress intensity factors computed for coarse mesh 


Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 
P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psi\/in 

Kn in psi'/in 

FEM 

634 

1 

1 

6.871998271 

-0.00186705222 

0.871918492 

-0.00213316604 

GFEM 

634 

1 

1 

0.764668694 

• -0.0017483637 

0.766995734 

-0.00157027579 

FEM 

634 

2 

1 

0.82840364 

-0.000819087014 

0.828413002 

-0.000812332551 

GFEM 

634 

2 

1 

0.978282491 

-0.00123810708 

0.978622215 

-0.00106178627 

Tabl 

e 5.2: Stress intensity factors computed for refined mesh 

Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi'/in 

Kn psi'/in 

FEM 

2564 

1 

1 

0.626637336 

-0.000337844132 

0.626767391 

-0.000283837087 

GFEM 

2564 

1 

1 

0.500493715 

-0.000460228968 

0.50052228 

-0.000217512687 

FEM 

2564 

1 

2 

0.849307746 

-0.000258345046 

0.849500731 

-0.000256760112 

GFEM 

2564 

1 

2 

0.855961206 

-0.00104771529 

0.855961206 

-0.000787968813 


Infinite Plate in Mode I Loading 


Ki = 0.895 


■psiy in 


Kn = 0 psiyin 


Discussion 

1. The theoretical value of Kj and Kn for finite plate, after applying finite plate corrections, 
are Ki = 0.895 and Kn = 0. The value obtained by GFEM and FEM with either 
p = 2 and one layer or p = 1 and two layer are sufficiently close to the theoretical 

values. 
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2. The contour integral method applied in the study is highly sensitive to the distance of the 
contour from the crack tip. As the contour is brought closer to the crack tip, the quality 
of the extracted stress intensity factor deteriorates. 

3. The stiffness intensity factors obtained at the left crack tip and right crack tip do not match 
exactly because of the lack of symmetry in the mesh used. Though, it should be noted 
that the discrepancy is insignificant. 

4. For coarse mesh, the results are not good beacause the crack length is too small. The 
contour integral employed in the extraction of the stress intensity factors will be severely 
aff’ected by the neighbouring crack tip. This is the reason that Fleming et al [4] suggested 
mesh refinement at the crack root for getting better results. 

Problem 2 : Centrally located straight edge long crack 

In this problem a centre cracked plate with the following parameters is analysed : 

Domain = 5 x 3.5 in^ 

Crack Length = 1.0 m 

A different crack length is chosen to demonstrate the effect of increasing crack length on the 
solution and how the higher order shape function affects the solution. In this problem the mesh 
employed is coarsh(see figure 5.5) The computed values of the stress intensity factors are reported 
in table 5.3 and 5.4 



Figure 5.5; Straight edge crack of length of lin 


Table 5.3: Stress intensity factors computed for p = 1 


Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi'/in 

Kti psi\/in 

FEM 

628 

1 

1 

0.90438288 

-0.00278504271 

0.904393975 

-0.0028196911 

GFEM 

628 

1 

1 

0.744402238 

-0.00184098385 

0.74466355 

-0.00191254745 

FEM 

628 

1 

2 

0.764668694 

-0.00236289666 

0.766995734 

-0.0023754772 

GFEM 

628 

1 

2 

1.25931867 

-0.00262408054 

1.26092989 

-0.00329535918 


Table 5.4: Stress intensity factors computed for p = 2 


Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki m psi'/in 

Kji psi'/in 

FEM 

628 

2 

1 

0.841305802 

-0.0017483637 

0.841326947 

-0.00157027579 



2 

1 

0.969028848 

-0.001050745056 

0.968735869 

-0.00149841222 

FEM 

628 

2 

2 

1.21598151 

-0.0017483637 

1.21636486 

-0.00157027579 



2 

2 

1.27639092 

-0.00104771529 

1.27570483 

-0.000787968813 


Infinite plate in Mode I loading 


Ki = 1.28 psivin 


Kii = 0 psi\/ in 


Discussion 


1. For p = 1, when the contour corresponds to two layer of elements, GFEM is very accurate 
while the standard FEM is not. Further using one layer of elements around the crack tip 
seems to be inefficient for the extraction of stress intensity factors. 

2. For p = 2, GFEM is clearly far superior to standard FEM. 
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3. Note that in this problem, the interference from the neighbouring crack tip is very small 
leading to more accurate extracted quantities. 

4. Inspite of a coarse unsymmetric mesh, the left and right stress intensity factor have in- 
significant mismatch. 


Problem 3 : Arbitrary shaped crack 

An arbitrary shaped crack with the following parameters is analysed : 


Domain 
Crack Segment : 

nodel 

node2 

nodes 

node4 

The plate with the crack is shown in figure 5.6. 1 
in table 5.5. 


= 5 X 3.5 in 

= (2.00,1.70) 

= (2.25,1.70) 

= (2.50,1.80) 

= (2.75,1.80) 

re computed stress intensity factors are reported 



Figure 5.6: Problem 3 


Table 5.5: Stress intensity factors computed for different p 


Scheme 

No. of 
. Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi\/in 

Kjj psi\/In 

GFEM 

636 

1 

1 

0.990401579 

0.0154856607 

0.982036791 

0.010429162 

GFEM 

636 

2 

1 

1.2372214 

0.0320909648 

1.24435412 

0.0309407329 

GFEM 

636 

3 

1 

1.20827447 

0.0536397577 

1.19315666 

0.022451595 


Discussion 

1. The stress intensity factors calculated for p = 2 and p = 3 are quite close and seems 
to converge (i.e. discrepancy is less than 5 pc). 

2. Note that the Kjj is essentially zero as compared to Kj and the crack will tend to grow 
straight. The left and the right stress intensity factor do not match because mesh is crude 
and unsymmetric, however, the discrepancy is minimal. 

3. To improve the result, one or two level of mesh refinemnet near the crack tip is desirable. 

Problem 4 : Arbitrary shaped crack 

An arbitrary shaped crack with the following parameters is analysed : 


Domain = 

5 X 3.5 

Segment : 

nodel = 

(2.00,1.70) 

node2 = 

(2.25,1.72) 

node3 = 

(2.50, 1.75) 

node4 = 

(2.75, 1.70) 

node5 = 

(3.00, 1.75) 


The arbitrary shaped crack is modelled as shown in figure 5.7. The computed stress intensity 
factors are reported in table 5.6 and 5.7. 
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Figure 5.7: Problem 4 


Table 5.6: Stress intensity factors computed for one layer 


Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P , 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Kj in psi\/in 

Kji psi\/in 

GFEM 

584 

1 

1 

1.13418414 

-0.00521811786 

1.08020361 

-0.0412967273 

GFEM 

584 

2 

1 

1.28536871 

0.0419222697 

1.28719329 

-0.063198877 

GFEM 

584 

3 



1 

1.47696526 

-0.0364505828 

1.41330589 

0.027973453 

Tal 

Die 5.7: Stress intensity fact 

,ors computed for two layer 

Scheme 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki in psWin 

Kii psiVin 

GFEM 

584 

1 

2 

1.5457882 

-0.051235026 

1.56144341 

0.0193207251 



2 

2 

1.4236723 

-0.0542374126 

1.38937786 

0.0665384271 
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Discussion 

1. The results for 1 layer of elements is quite poor. However for p = 3 the results seem to 
be close to that obtained from extraction using two layer of elements around the crack tip. 

2. Again in this case Ku is negligible because the crack behaves like a centrally placed crack. 


Problem 5 : Crack growth for two symmetrically located crack 

The growth of two symmetrically placed straight cracks with the following initial parameters is 
analysed here and the corresponding position of the cracks are shown in figure (5.8 - 5.11). 


Domain = 
Crack Profile 1 : 

nodel = 
node2 = 
nodes = 
node4 = 
node5 = 
Crack Profile 2 : 

nodel = 
node2 = 
nodes = 
node4 = 
node5 = 


5 X S.5 in^ 

(1.500. 1.75) 

(1.625.1.75) 

(1.750. 1.75) 

(1.875.1.75) 
(2.000,1.75) 

(S.OOO, 1.75) 

(5.125.1.75) 

(5.250. 1.75) 

(5.575. 1.75) 

(5.500. 1.75) 


The purpose of the present problem is to estimate the robustness of the code for moving geom- 
etry. Here in each step, the crack length gets increased and the analysis is repeated. The crack 
extension direction is determined according to maximum tangential stress criterion detailed ear- 
lier in chapter 2. The extracted stress intensity factors for different p values are reported in table 
5.8-5.1S. 


36 



step 1 



Figure 5.8: Problem 5 crack growth step 1 


Table 5.8: S1 

tress intensity factors computed at the end of ste] 

1 with p = 1 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Kj in psiVin 

Kji psi'/in 

1 

2564 

1 

2 

0.827900965 

-0.000352688633 

0.85007153 

-0.000317667391 

2 


1 

2 

0.850005244 

-0.000285569277 

0.827116411 

-0.0000478630334 

Table 5.9: S1 

:ress intensity factors computed at the end of ste] 

1 with p = 2 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Ki in psiVin 

Kii psi\/in 


2562 

2 

2 

0.853412354 

-0.000240084291 

0.869690615 

-0.000255675898 

2 

2562 

2 

2 

0.869697323 

-0.000280393539 

0.853348802 

-0.0000264159494 
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Step 2 



Figure 5.9: Problem 5 crack growth step 2 


Table 5.10: Stress intensity factors computed at the end of step 2 with p = 1 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiy/in 

Kji psi\/in 

1 

2564 

1 

2 

0.820305226 

-0.000301340611 



2 

2564 

1 

2 

0.865398794 

-0.000570112828 

0.819049369 

0.00263806512 


Table 5.11: Stress intensity factors computed at the end of step 2 with p = 2 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Kj in psiy/i^ 

Kji psiy/in 

1 

2560 

2 

2 

0.848856914 

-0.000126843033 

0.883541027 

-0.000481974206 

2 


2 

2 

0.881666649 

-0.000156524428 

0.85113107 

-0.000796571682 
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step 3 



Figure 5.10: Problem 5 crack growth step 3 


Table 5.12: Stress intensity factors computed at the end of step 3 with p = 1 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiVin 

Kir psiy/in 

1 

2510 

1 

i 

2 

0.912694996 

-0.000292969031 

0.964575631 

-0.00703674212 

2 

2510 

1 

2 

0.918500545 

-0.00261285263 

0.829219297 

-0.00526810027 


Table 5.13: Stress intensity factors computed at the end of step 3 with p = 2 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Elements 
to compute 
Interaction 
Integral 

Kr in psiy/in 

Kii psiyfin 

1 

2566 

2 

2 

0.911692321 

-0.000531385793 

1.00692173 

-0.00147615579 


2566 

2 

2 

1.00926826 

-0.000218973043 

0.921226773 

0.00378581801 
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Final step 



Figure 5.11: Problem 5 crack growth step 4 


Discussion 

1. Following the theoretical analysis, geometry and location of the crack, both the crack should 
propagate along a straight line. This is what we get out of our numerical analysis. 

2. As the crack grows, the interaction between two cracks become more severe leading to an 
increase in stress intensity factor. 

3. The crack propagation is stopped prior to merger because of the limitations of the current 
meshing process. 

Problem 6 : crack growth for two airbitrarily located crack 

In this problem two cracks are considered. However the cracks are located off center as shown 
in the figure(5.12). The details of the domain and the crack geometry are as follows : 

Domain =5 x 3.5 in^ 

Crack Profile 1 : 
nodel = (1.500, 1.25) 
node2 = (1.625, 1.25) 
nodes = (1.750,1.25) 
node4 = (1.875, 1.25) 
nodes = (2.000, 1.25) 

Crack Profile 2 : 
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nodel = (3.000,2.25) 
node2 = (3.125,2.25) 
node3 = (3.250,2.25) 
node4 = (3.375,2.25) 
node5 = (3.500, 2.25) 

Each of the crack segment is taken as a series of line segments joining the above 5 nodes. After 
each step of iteration, the crack profiles are given proper crack extension in the right direction 
and the whole analysis is repeated. The computed stress intensity factors at the end of each step 
are reported in the table 5.14 - 5.23. The corresponding crack geometries at the beginning of 
each step are shown in figure (5.12 - 5.21). 

Step 1 
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Figure 5.12: Problem 5 crack growth step 1 


Table 5.14: Stress ini 

tensity factors computed at the end of step 1 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psiVin 

Kn psiV in 

1 

2584 

1 

2 

0.835341791 

0.0121879137 

0.850736133 

0.0113071723 

2 

2584 

1 

2 

0.850696036 

0.0113071723 


0.0121700265 
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step 2 



Figure 5.13: Problem 6 crack growth step 2 


Table 5.15: Stress inl 

tensity factors computed at the end of step 2 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiVin 

Kn psiVin 

1 

2596 

1 

2 

0.830585879 

0.00567496787 

0.851276091 

-0.00508734123 


2596 

t 

1 

2 

0.864274472 

0.00791880615 

0.828841718 

0.00307971677 
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step 3 



Figure 5.14: Problem 6 orack growth step 3 


Table 5.16: Stress intensity factors computed at the end of step 3 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiy/in 

Kri psiVin 

1 

2578 

1 

2 

1.00135083 

-0.00158095561 

1.09011459 

-0.0132891375 

2 

2578 

1 

2 

1.093026525 

0.00643484579 

1.02283768 

0.00345989243 
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step 4 
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Figure 5.15: Problem 6 crack growth step 4 


Table 5.17: Stress in1 

:ensity factors computed at the end of step 4 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psi\fin 

Kji psi'/in 

1 

2570 

1 

2 

1.36818268 

-0.0211871321 

1.48479463 

-0.0986211241 

2 

2570 

1 

2 

1.4507246 

-0.0194111176 

1.36069751 

-0.000822044013 
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Figure 5.16: Problem 6 crack growth step 5 


Table 5.18: Stress intensity factors computed at the end of step 5 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psiy/^ 

JT/j psiy/in 

1 

2570 

1 

2 

1.71884708 

-0.0403886554 

1.86253646 

-0.0604739141 

to 


1 

2 

1.84575045 

-0.200539113 

1.66395861 

-0.117210394 
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step 6 



Figure 5.17: Problem 6 crack growth step 6 


Table 5.19: Stress ini 

tensity factors computed at the end of step 6 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiV^ 

Kn psiy/in 

1 

2582 

1 

'2 

2.02041082 

-0.0529359263 

2.19797146 

-0.313187923 

2 

2582 

1 

2 

2.126893 

-0.336561875 

2.09673951 

-0.0556087537 
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step 7 



Figure 5.18: Problem 6 crack growth step 7 


Table 5.20: Stress im 

tensity factors computed at the end of step 7 

Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psiy/in 

Kn psiVin 

1 

2580 

1 

2 

2.48713937 

-0.234273036 

2.56750505 

-0.308852834 

2 

2580 

1 

2 

2.43620877 

-0.309309869 

2.44884128 

-0.23661355 
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Figure 5.19: Problem 6 crack growth step 8 


Table 5.21: Stress intensity factors computed at the end of step 8 

lo. of 

Ele- 

aents 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psi\/in 

Kji psiy/in 




2.94236899 

-0.324194315 

i92 

1 

2 

2.6786357 

-0.315583614 




2.72784844 

-0.491364253 

i92 

1 

•2 

2.90369245 

-0.276572426 




step 9 



Figure 5.20: Problem 6 crack growth step 9 


Table 5.22: Stress intensity factors computed at the end of step 9 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Kj in psiV^ 

Kji psi\/in 

1 

2572 

1 

2 ^ 

3.30572489 

-0.402490819 

2.97397514 

-0.428100794 

2 

2572 

1 

2 

2.98920629 

-0.354454145 

3.22757431 

-0.445701305 
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step 10 



Figure 5.21: Problem 6 crack growth step 10 


Table 5.23: Stress intensity factors computed at the end of step 10 


Crack 

No. of 
Ele- 
ments 

Order of 
Shape 
Functions , 

P 

Layer of 
Element to 
compute 
Interaction 
Integral 

Ki in psis/in 

Kjj psiVin 

1 

2560 

1 

2 

3.65941563 

-0.513562382 

3.21282555 

-0.270626624 

2 

2560 

1 

2 

3.17263255 

-0.394161767 

3.71137782 

-0.449254288 


Discussion 

1. A close look at the problem reveals that the crack positions are unaltered by a rotation of 
the domain by 180°. So one would expect the cracks to propagate in a self similar manner. 
The numerical results demonstrate that indeed it is true. 
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2. As the cracks grow, the interaction between the two cracks increases leading to an increase 
in the observed stress intensity factors. 

3. The cracks tend to turn towards each other and towards the fixed boundary points at the 
two ends. This confirms with the observed photoelastic data for similar experiment as well 
as to the intuitive reasoning that near the boundary, the crack will tend to grow towards 
the region of high stress concentration. 

4. The mixed mode stress intensity factor is mainly Ki dominated. Mode II is important 
from the point of view of the orientation of the crack growth direction. The crack turns 
primarily due to Kji. 
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Chapter 6 
Closure 


6.1 Summary 

In the present work, we developed a general computational code to model crack growth. The 
following conclusions can be drawn based on the present work : 

1. The presence of multiple singularities in the domain has been successfully addressed. How- 
ever more detailed analysis is to be done to improve the results further. 

2. The generalized finite element is superior than the standard finite element solution. 

3. The present approach has been extended to arbitrary shaped crack. 

4. Mesh refinement is necessary in some cases, as in problem 1, where due to small crack 
length the local disturbance from the nearby cracktip is affecting the solution. 

5. The stress intensity factors calculated from the interaction integral computed along the 
two layer of elements around the crack tip yield better results. 

6. Higher order shape functions along with little mesh refinement could be very effective tool 
to characterise the crack parameters. 


6.2 Suggestions for future work 

1. The major difficulty in crack growth modelling in any of the finite element approach is 
to generate a mesh at each step conforming to the crack surface. The present study 
uses automatic mesh generator in crude form. The mesh generator should be improved by 
incorporating sub-domain wise meshing. The whole domain could be divided into a number 
of sub domains such that the crack surface falls on one of the boundary. This excludes the 
possiblity of any internal cut-outs. 

gpnfcg 9^0 A 
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2. In the present work, we have enriched only the node sitting at the crack tip. This could 
be extended to enrich the whole zone around the crack tip. This will improve the results 
immensely. 

3. The enrichment functions used in the present method is as proposed by Fleming et al [4]. 
However, much better results could be obtained by using the exact functions as derived for 
re-entrant corner. 

4. Mesh refinement might be necessary in some cases particularly when the cracks come 
closer to each other. This would improve the results immensely. Instead of mesh refinement 
throughout the whole domain, it is better to refine only those elements which are associated 
with the node sitting on the crack tip. 

5. In the present study contour form of interaction integral is applied to extract stress intensity 
factor in mixed mode analysis. However, domain form of interaction integral based on cut- 
off function method (see [12]) or contour form of interaction integral around an annular 
ring surrounding the crack tip (see [12]) would yield better results. 
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